startdatenum = 1817;
enddatenum = 1914;

data = readtable('../../data/Netherlands/NetherlandsPF_1815_1938.xlsx');
date = data.year;

startdate = find(date == startdatenum);
enddate = find(date == enddatenum);
T = length(date(startdate:enddate));

option.detrend = 1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Load data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

time = data.Var1(startdate:enddate);
pi = data.pi(startdate:enddate);
rpcgdpgr = data.x(startdate:enddate);
tau = data.tau(startdate:enddate);
g = data.g(startdate:enddate);
ynom1 = log(1 + data.y1(startdate:enddate));
ynom10 = log(1 + data.y10(startdate:enddate));
surplus = tau - g;

debt = data.debt_gdp(startdate:enddate) / 100;

deltalogg = [log(data.g(startdate)) - log(data.g(startdate - 1)); log(g(2:end)) - log(g(1:end - 1))];
deltalogtau = [log(data.tau(startdate)) - log(data.tau(startdate - 1)); log(tau(2:end)) - log(tau(1:end - 1))];

yieldspr = ynom10 - ynom1;

x0 = mean(rpcgdpgr);
pi0 = mean(pi);
tau0 = mean(deltalogtau);
g0 = mean(deltalogg);
y0nom_1 = mean(ynom1);
y0nom_10 = mean(ynom10);
surplus0 = mean(surplus);
yspr0 = mean(yieldspr);

%%%debt market value adjustment based on van Riel (2018) Chapter 7
slope =- (1/3 - 760.1 / (760.1 + 837.0)) / (1828 - 1815);
ratio = 1/3 + (0:1:97)' * slope - 2 * slope;
ratio(ratio > 1) = 1;
debt_adjusted = (debt .* ratio) ./ ynom10 * 0.025 + debt .* (1 - ratio) * 0.019;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Detrend tax and spending
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if option.detrend == 1
    tts(1) = tau(1);
    gts(1) = g(1);

    for t = 2:T
        gts(t) = gts(t - 1) * exp(deltalogg(t) - g0);
        tts(t) = tts(t - 1) * exp(deltalogtau(t) - tau0);
    end

end

taxrevgdp = tau;
spendgdp = g;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Assembly
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% define the ordering in the VAR
inflpos = 1;
y1pos = 2;
ysprpos = 3;
gdppos = 4;
dtpos = 5;
coint_tax = 6;
dgpos = 7;
coint_spending = 8;

%% raw data
clear Xraw;
Xraw(:, inflpos) = pi;
Xraw(:, gdppos) = rpcgdpgr;
Xraw(:, y1pos) = ynom1;
Xraw(:, ysprpos) = yieldspr;
Xraw(:, dtpos) = deltalogtau;
Xraw(:, coint_tax) = log(tau);
Xraw(:, dgpos) = deltalogg;
Xraw(:, coint_spending) = log(g);

% VAR elements
infl = pi - pi0;
x = rpcgdpgr - x0;
yield1 = ynom1 - y0nom_1;
nomyspr = yieldspr - yspr0;
dg = deltalogg - g0;
dt = deltalogtau - tau0;

X2(:, inflpos) = infl;
X2(:, y1pos) = yield1;
X2(:, ysprpos) = nomyspr;
X2(:, gdppos) = x;
X2(:, dtpos) = dt;
X2(:, coint_tax) = log(t) - mean(log(t));
X2(:, dgpos) = dg;
X2(:, coint_spending) = log(g) - mean(log(g));

% trend adjustment
if option.detrend == 1
    X2(:, coint_tax) = log(tts) - mean(log(tts));
    X2(:, coint_spending) = log(gts) - mean(log(gts));
end

z_var = X2(2:end, :);
z_varlag = X2(1:end - 1, :);
Y = z_var;
X = z_varlag;

N = cols(X2);
I = eye(N);
I_pi = I(:, inflpos);
I_gdp = I(:, gdppos);
I_y1 = I(:, y1pos);
I_yspr = I(:, ysprpos);
I_dt = I(:, dtpos);
I_dg = I(:, dgpos);
I_coint_tax = I(:, coint_tax);
I_coint_spending = I(:, coint_spending);

% actual, non-detrended series of tax and spending
X2t = X2;
X2t(:, coint_tax) = log(tau) - mean(log(tau));
X2t(:, coint_spending) = log(g) - mean(log(g));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Estimate Psi and Sig
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Psi = zeros(N, N);
Psi_se = zeros(N, N);

for i = [inflpos y1pos ysprpos gdppos dtpos dgpos]
    regr = ols(Y(:, i), [ones(T - 1, 1), X(:, :)]);
    c(i) = regr.beta(1);
    c_se(i) = regr.bstd(1);
    Psi(i, :) = regr.beta(2:end)';
    Psi_se(i, :) = regr.bstd(2:end)';
    R2(i) = regr.rsqr;
    R2bar(i) = regr.rbar;
    eps(:, i) = Y(:, i) - c(i) - X * Psi(i, :)';
    clear regr;
end

Psi(coint_tax, :) = Psi(dtpos, :);
Psi(coint_spending, :) = Psi(dgpos, :);
Psi(coint_tax, coint_tax) = Psi(coint_tax, coint_tax) + 1;
Psi(coint_spending, coint_spending) = Psi(coint_spending, coint_spending) + 1;

Psi_se(coint_tax, :) = Psi_se(dtpos, :);
Psi_se(coint_spending, :) = Psi_se(dgpos, :);

Sigma = cov(eps(:, [inflpos y1pos ysprpos gdppos dtpos dgpos]));
tmp = chol(Sigma, 'lower');
Sig = zeros(N, N);
Sig([inflpos y1pos ysprpos gdppos], [inflpos y1pos ysprpos gdppos]) = ...
    tmp(1:4, 1:4);

Sig(dtpos, [inflpos y1pos ysprpos gdppos dtpos]) = tmp(5, 1:5);
Sig(dgpos, [inflpos y1pos ysprpos gdppos dtpos dgpos]) = tmp(6, 1:6);

Sig(coint_tax, :) = Sig(dtpos, :);
Sig(coint_spending, :) = Sig(dgpos, :);

eps = [eps(:, [inflpos y1pos ysprpos gdppos]), ...
                eps(:, dtpos), eps(:, dtpos), ...
                eps(:, dgpos), eps(:, dgpos)];
Sigma = Sig * Sig';

tstat = Psi ./ Psi_se;
